// ---------------------------------------------------------------------
//
// Copyright (c) 2014 - 2019 by the IBAMR developers
// All rights reserved.
//
// This file is part of IBAMR.
//
// IBAMR is free software and is distributed under the 3-clause BSD
// license. The full text of the license can be found in the file
// COPYRIGHT at the top level directory of IBAMR.
//
// ---------------------------------------------------------------------

/////////////////////////////// INCLUDES /////////////////////////////////////

#include "ibtk/ibtk_utilities.h"
#include "ibtk/muParserRobinBcCoefs.h"
#include "ibtk/namespaces.h" // IWYU pragma: keep

#include "ArrayData.h"
#include "BoundaryBox.h"
#include "Box.h"
#include "CartesianGridGeometry.h"
#include "CartesianPatchGeometry.h"
#include "Index.h"
#include "IntVector.h"
#include "Patch.h"
#include "tbox/Array.h"
#include "tbox/Database.h"
#include "tbox/Pointer.h"
#include "tbox/Utilities.h"

#include "muParser.h"
#include "muParserError.h"

#include <algorithm>
#include <array>
#include <map>
#include <ostream>
#include <string>
#include <utility>
#include <vector>

namespace SAMRAI
{
namespace hier
{
template <int DIM>
class Variable;
} // namespace hier
} // namespace SAMRAI

/////////////////////////////// NAMESPACE ////////////////////////////////////

namespace IBTK
{
/////////////////////////////// STATIC ///////////////////////////////////////

namespace
{
static const int EXTENSIONS_FILLABLE = 128;
}

/////////////////////////////// PUBLIC ///////////////////////////////////////

muParserRobinBcCoefs::muParserRobinBcCoefs(const std::string& object_name,
                                           Pointer<Database> input_db,
                                           Pointer<CartesianGridGeometry<NDIM> > grid_geom)
    : d_grid_geom(grid_geom)
{
#if !defined(NDEBUG)
    TBOX_ASSERT(!object_name.empty());
    TBOX_ASSERT(input_db);
#else
    NULL_USE(object_name);
#endif
    // Read in user-provided constants.
    Array<std::string> db_key_names = input_db->getAllKeys();
    for (int k = 0; k < db_key_names.size(); ++k)
    {
        const std::string& name = db_key_names[k];
        if (input_db->isDouble(name))
        {
            d_constants[name] = input_db->getDouble(name);
        }
        else if (input_db->isFloat(name))
        {
            d_constants[name] = input_db->getFloat(name);
        }
        else if (input_db->isInteger(name))
        {
            d_constants[name] = input_db->getInteger(name);
        }
    }

    // Initialize the parsers with data read in from the input database.
    for (int d = 0; d < 2 * NDIM; ++d)
    {
        std::string key_name;
        const std::string postfix = "_function_" + std::to_string(d);

        key_name = "acoef" + postfix;
        if (input_db->isString(key_name))
        {
            d_acoef_function_strings.push_back(input_db->getString(key_name));
        }
        else
        {
            d_acoef_function_strings.push_back("0.0");
            TBOX_WARNING("muParserRobinBcCoefs::muParserRobinBcCoefs():\n"
                         << "  no function corresponding to key ``" << key_name << "'' found for side = " << d
                         << "; using acoef = 0.0." << std::endl);
        }
        try
        {
            d_acoef_parsers[d].SetExpr(d_acoef_function_strings.back());
        }
        catch (mu::ParserError& e)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  error: " << e.GetMsg() << "\n"
                       << "  in:    " << e.GetExpr() << "\n");
        }
        catch (...)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  unrecognized exception generated by muParser library.\n");
        }

        key_name = "bcoef" + postfix;
        if (input_db->isString(key_name))
        {
            d_bcoef_function_strings.push_back(input_db->getString(key_name));
        }
        else
        {
            d_bcoef_function_strings.push_back("0.0");
            TBOX_WARNING("muParserRobinBcCoefs::muParserRobinBcCoefs():\n"
                         << "  no function corresponding to key ``" << key_name << "'' found for side = " << d
                         << "; using bcoef = 0.0." << std::endl);
        }
        try
        {
            d_bcoef_parsers[d].SetExpr(d_bcoef_function_strings.back());
        }
        catch (mu::ParserError& e)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  error: " << e.GetMsg() << "\n"
                       << "  in:    " << e.GetExpr() << "\n");
        }
        catch (...)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  unrecognized exception generated by muParser library.\n");
        }

        key_name = "gcoef" + postfix;
        if (input_db->isString(key_name))
        {
            d_gcoef_function_strings.push_back(input_db->getString(key_name));
        }
        else
        {
            d_gcoef_function_strings.push_back("0.0");
            TBOX_WARNING("muParserRobinBcCoefs::muParserRobinBcCoefs():\n"
                         << "  no function corresponding to key ``" << key_name << "'' found for side = " << d
                         << "; using gcoef = 0.0." << std::endl);
        }
        try
        {
            d_gcoef_parsers[d].SetExpr(d_gcoef_function_strings.back());
        }
        catch (mu::ParserError& e)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  error: " << e.GetMsg() << "\n"
                       << "  in:    " << e.GetExpr() << "\n");
        }
        catch (...)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  unrecognized exception generated by muParser library.\n");
        }
    }

    // Define the default and user-provided constants.
    std::vector<mu::Parser*> all_parsers(3 * 2 * NDIM);
    for (int d = 0; d < 2 * NDIM; ++d)
    {
        all_parsers[3 * d] = &d_acoef_parsers[d];
        all_parsers[3 * d + 1] = &d_bcoef_parsers[d];
        all_parsers[3 * d + 2] = &d_gcoef_parsers[d];
    }
    const double pi = 3.1415926535897932384626433832795;
    const double* const xLower = grid_geom->getXLower();
    const double* const xUpper = grid_geom->getXUpper();
    for (const auto& parser : all_parsers)
    {
        // Various names for pi.
        parser->DefineConst("pi", pi);
        parser->DefineConst("Pi", pi);
        parser->DefineConst("PI", pi);

        // The extents of the domain.
        for (unsigned int d = 0; d < NDIM; ++d)
        {
            const std::string postfix = std::to_string(d);

            parser->DefineConst("X_LOWER" + postfix, xLower[d]);
            parser->DefineConst("X_lower" + postfix, xLower[d]);
            parser->DefineConst("x_lower" + postfix, xLower[d]);
            parser->DefineConst("x_LOWER" + postfix, xLower[d]);
            parser->DefineConst("X_Lower" + postfix, xLower[d]);
            parser->DefineConst("X_lower" + postfix, xLower[d]);
            parser->DefineConst("XLower" + postfix, xLower[d]);
            parser->DefineConst("Xlower" + postfix, xLower[d]);
            parser->DefineConst("x_Lower" + postfix, xLower[d]);
            parser->DefineConst("x_lower" + postfix, xLower[d]);
            parser->DefineConst("xLower" + postfix, xLower[d]);
            parser->DefineConst("xlower" + postfix, xLower[d]);

            parser->DefineConst("X_LOWER_" + postfix, xLower[d]);
            parser->DefineConst("X_lower_" + postfix, xLower[d]);
            parser->DefineConst("x_lower_" + postfix, xLower[d]);
            parser->DefineConst("x_LOWER_" + postfix, xLower[d]);
            parser->DefineConst("X_Lower_" + postfix, xLower[d]);
            parser->DefineConst("X_lower_" + postfix, xLower[d]);
            parser->DefineConst("XLower_" + postfix, xLower[d]);
            parser->DefineConst("Xlower_" + postfix, xLower[d]);
            parser->DefineConst("x_Lower_" + postfix, xLower[d]);
            parser->DefineConst("x_lower_" + postfix, xLower[d]);
            parser->DefineConst("xLower_" + postfix, xLower[d]);
            parser->DefineConst("xlower_" + postfix, xLower[d]);

            parser->DefineConst("X_UPPER" + postfix, xUpper[d]);
            parser->DefineConst("X_upper" + postfix, xUpper[d]);
            parser->DefineConst("x_upper" + postfix, xUpper[d]);
            parser->DefineConst("x_UPPER" + postfix, xUpper[d]);
            parser->DefineConst("X_Upper" + postfix, xUpper[d]);
            parser->DefineConst("X_upper" + postfix, xUpper[d]);
            parser->DefineConst("XUpper" + postfix, xUpper[d]);
            parser->DefineConst("Xupper" + postfix, xUpper[d]);
            parser->DefineConst("x_Upper" + postfix, xUpper[d]);
            parser->DefineConst("x_upper" + postfix, xUpper[d]);
            parser->DefineConst("xUpper" + postfix, xUpper[d]);
            parser->DefineConst("xupper" + postfix, xUpper[d]);

            parser->DefineConst("X_UPPER_" + postfix, xUpper[d]);
            parser->DefineConst("X_upper_" + postfix, xUpper[d]);
            parser->DefineConst("x_upper_" + postfix, xUpper[d]);
            parser->DefineConst("x_UPPER_" + postfix, xUpper[d]);
            parser->DefineConst("X_Upper_" + postfix, xUpper[d]);
            parser->DefineConst("X_upper_" + postfix, xUpper[d]);
            parser->DefineConst("XUpper_" + postfix, xUpper[d]);
            parser->DefineConst("Xupper_" + postfix, xUpper[d]);
            parser->DefineConst("x_Upper_" + postfix, xUpper[d]);
            parser->DefineConst("x_upper_" + postfix, xUpper[d]);
            parser->DefineConst("xUpper_" + postfix, xUpper[d]);
            parser->DefineConst("xupper_" + postfix, xUpper[d]);
        }

        // User-provided constants.
        for (const auto& constant : d_constants)
        {
            parser->DefineConst(constant.first, constant.second);
        }

        // Variables.
        parser->DefineVar("T", &d_parser_time);
        parser->DefineVar("t", &d_parser_time);
        for (unsigned int d = 0; d < NDIM; ++d)
        {
            const std::string postfix = std::to_string(d);
            parser->DefineVar("X" + postfix, d_parser_posn.data() + d);
            parser->DefineVar("x" + postfix, d_parser_posn.data() + d);
            parser->DefineVar("X_" + postfix, d_parser_posn.data() + d);
            parser->DefineVar("x_" + postfix, d_parser_posn.data() + d);
        }
    }
    return;
} // muParserRobinBcCoefs

void
muParserRobinBcCoefs::setBcCoefs(Pointer<ArrayData<NDIM, double> >& acoef_data,
                                 Pointer<ArrayData<NDIM, double> >& bcoef_data,
                                 Pointer<ArrayData<NDIM, double> >& gcoef_data,
                                 const Pointer<Variable<NDIM> >& /*variable*/,
                                 const Patch<NDIM>& patch,
                                 const BoundaryBox<NDIM>& bdry_box,
                                 double fill_time) const
{
    const Box<NDIM>& patch_box = patch.getBox();
    const hier::Index<NDIM>& patch_lower = patch_box.lower();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom = patch.getPatchGeometry();

    const double* const x_lower = pgeom->getXLower();
    const double* const dx = pgeom->getDx();

    // Loop over the boundary box and set the coefficients.
    const unsigned int location_index = bdry_box.getLocationIndex();
    const unsigned int bdry_normal_axis = location_index / 2;
    const Box<NDIM>& bc_coef_box =
        (acoef_data ? acoef_data->getBox() :
                      bcoef_data ? bcoef_data->getBox() : gcoef_data ? gcoef_data->getBox() : Box<NDIM>());
#if !defined(NDEBUG)
    TBOX_ASSERT(!acoef_data || bc_coef_box == acoef_data->getBox());
    TBOX_ASSERT(!bcoef_data || bc_coef_box == bcoef_data->getBox());
    TBOX_ASSERT(!gcoef_data || bc_coef_box == gcoef_data->getBox());
#endif

    const mu::Parser& acoef_parser = d_acoef_parsers[location_index];
    const mu::Parser& bcoef_parser = d_bcoef_parsers[location_index];
    const mu::Parser& gcoef_parser = d_gcoef_parsers[location_index];
    d_parser_time = fill_time;
    for (Box<NDIM>::Iterator b(bc_coef_box); b; b++)
    {
        const hier::Index<NDIM>& i = b();
        for (unsigned int d = 0; d < NDIM; ++d)
        {
            if (d != bdry_normal_axis)
            {
                d_parser_posn[d] = x_lower[d] + dx[d] * (static_cast<double>(i(d) - patch_lower(d)) + 0.5);
            }
            else
            {
                d_parser_posn[d] = x_lower[d] + dx[d] * (static_cast<double>(i(d) - patch_lower(d)));
            }
        }
        try
        {
            if (acoef_data) (*acoef_data)(i, 0) = acoef_parser.Eval();
            if (bcoef_data) (*bcoef_data)(i, 0) = bcoef_parser.Eval();
            if (gcoef_data) (*gcoef_data)(i, 0) = gcoef_parser.Eval();
        }
        catch (mu::ParserError& e)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  error: " << e.GetMsg() << "\n"
                       << "  in:    " << e.GetExpr() << "\n");
        }
        catch (...)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  unrecognized exception generated by muParser library.\n");
        }
    }
    return;
} // setBcCoefs

IntVector<NDIM>
muParserRobinBcCoefs::numberOfExtensionsFillable() const
{
    return EXTENSIONS_FILLABLE;
} // numberOfExtensionsFillable

/////////////////////////////// PROTECTED ////////////////////////////////////

/////////////////////////////// PRIVATE //////////////////////////////////////

/////////////////////////////// NAMESPACE ////////////////////////////////////

} // namespace IBTK

//////////////////////////////////////////////////////////////////////////////
